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ABSTRACT 

Most detected planet-bearing binaries are in wide orbits, for which a high inclination, ig, between 
the binary orbital plane and the plane of the planetary disk around the primary is likely to be 
common. In this paper, we investigate the intermediate stages - from planetcsimals to planetary 
embryos/cores - of planet formation in such highly inclined cases. Our focus is on the effects of 
gas drag on the planetesimals' orbital evolution, in particular on the evolution of the planetesimals' 
semimajor axis distribution and their mutual relative velocities. We first demonstrate that a non- 
evolving axisymmetric disk model is a good approximation for studying the effects of gas drag on 
a planetesimal in the highly inclined case (30° < is < 150°). We then find that gas drag plays a 
crucial role, and the results can be generally divided into two categories, i.e., the Kozai-on regime 
and the Kozai-off regime, depending on the specific value of %b- For both regimes, a robust outcome 
over a wide range of parameters is that, planetesimals migrate/jump inwards and pile up, leading to a 
severely truncated and dense planetesimal disk around the primary. In this compact and dense disk, 
collision rates are high but relative velocities arc low, providing conditions which are favorable for 
planetesimal growth, and potentially allow for the subsequent formation of planets. 
Subject headings: methods: numerical — planets and satellites: formation 



1. INTRODUCTION 

Planets in binaries could be common as stars 
are usually born in binary or multi-stars sys- 
tems (iDuauennov fc Mavorl 11991b iMathieu etHI 120001 : 
iDuchene et all I2007T ). Although current observations 
show that the multiplicit y-rate of the detected exoplane t 
host stars is around 17% (|Mugrauer fc Neuhauserll2009D . 
this fraction should be considered as a lower limit because 
of noticea ble selection effects against binaries in planet 
searches (Eggcnbe rgerl I2010T ). Most of the currently 
known planet-beari n g bin ar y systems are wide S -types 
(lEggenberger et all 120041 [Pesidera fc Barbieril 120071: 
lHaghighipourl 120101 ) . meaning the companion star acts 
as a distant satellite, typically orbiting the inner star- 
planet system over 100 AU away. Nevertheless, there are 
currently four systems with smaller separation of ~ 20 
AU, including the 7 Cephei (lHatzes et all 12003ft. G J 86, 
(jQueloz et al.ll2000D. H P 41004 dZucker et al.ll200l and 



constraints for theories of planet formation. Cur- 
rent theories of the core-accretion model of planet 
formation in single-planet systems (jLissauerl 119931 : 
lArmitagd l20Tot) envisage formation proceeding via a 
number of successive stages. In the first step, sub- 
mm dust grains are converted into 0.1-1000 km- 
sized planetesimals, the specific mechanism for which 
may involve (i) some form of grain coagulation 
(iWeidenschilling fc Cuzzil 119931: IWeidenschillingl 11997b 



Uohansen et al.l 12008b iTeiser fc Wurml 120091 ) and /or 
(ii) some form of fragmentation due to instabili- 
ties in the solid sub-com p onent of the protoplan- 
etary gas disk (iSafronovl 11969b iGoldreich fc Ward! 
19731: lYoudin fc Shul 120021: lYoudin fc Goodman! 12005b 
Youdin fc Johansenl 12007b iCuzzi et al.l 120081) . Reviews 



HP 196885 (ICorreia et al.l[2008tlChauvin et a.1 .11 2 1 OT) . In 
addition to the planets in circumprimary (S-type) orbits 
discussed above, planets in circu mbinary (P-type) orbits 
have been foun d in two systems (jSigurdsson et al.l 12003b 
iLee et aDl2009f ) . While it is possible that the recent de- 
tections of numero us short period eclipsing binaries by 
t he Kepler mission (jPrsa et alJlMTblSlawson et al.|[20Tll) 
may promote additional such circumbinary detections, 
their short periods are likely to preclude their use in the 
characterization of the circumprimary planets on which 
we concentrate in this investigation. 

The properties and diversities of planet-bearing 
double-star systems provide additional challenges and 

xicjiwci@gmail.com 



of these topics ar e given by iBlum fc Wurmt ([2008D and 
IChiang fc Youdinl ([2010D . In the second stage, mutual 
(planetesimal-planetesimal) collisions commence and the 
planetesimals grow yia runaway and oligarchic phases 
(IWetherill fc Stewart] Il989l: IKokubo fc Idal [19961 [1991 
lRafikovl2*0 03. 2004) to form protoplanetary embryos. Pe- 
pending on the birth-rate and birth size of planetesimals, 
there could be a transition phase— " snowball" phase— be- 
tween the first and second stages (jXie et al.ll2010bl ). In 
the third stage, giant impacts between these embryos 
commence, leading to the form ation of full-size terres- 
trial planets and gas-giant cores ([Chambers fe Wetherilll 
[T9^lLCTds"on fc Agnorll2"00l IKokubo et al.ll2006fl ~ 

Perturbative effects from stellar companions render the 
planet formation process in binary systems even more 
complex than that in single-star systems. For the first 
stage of planet formation outlined above, it is still un- 
clear whether and/or how the formation of planetesimals 
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from dust is affected by the binary star ([Pascucci et al.l 
120081 iZsom et al.l 120101 ): whereas for the third stage, 
it is thought that plane t formation is probably unaf- 
fected or even hastened (iBarbieri et al.l 120021: iQuintanal 
20041: iQuintana fc Lissauerl 120061 : iQuintana et a ' .1 120071 : 
Guedes et al.1 120081 ) in a region that extends to the 
outer limit for orbital stability, i.e., out to ^3-5 AU 
for a typical tight binar y of semi-major axis 20 AU 
(IHolman fc Wiegertl[l999l) . 

However, the intermediate (second) stage of the forma- 
tion process is much more problematic, as the outcome 
of planetesimal-planetesimal collisions is highly sensitive 
to the rela tive velocities between t he bodies during any 
encou nter ([Benz fc Asphaud 119991 : iStewart fc Leinhardtl 
l2009t) . The perturbations from a binary companion can 
excite planctesimal orbits and increase their mutual im- 
pact velocity, AV, to values that might exceed their es- 
cape velocity or ev en the critical velocity for the onset of 
erodi ng collisions ([Heppenheimerl 119781 : iWhitmire et al.l 
[19981) . 

Some earlier investigations ([Marzari fc Scholll l2000h 
suggested that the combination of binary perturbations 
and gas damping could force a strong pericenter align- 
ment between p lanetesimal orbi t s, red ucing AV^. How- 
ever the work of lThebault et all (120061 ) showed that it is 
only similarly sized objects that experience this reduc- 
tion, as in general the preferred pericenter alignment is 
strongly size-dependent, leading to an increasing AV for 
objects of different sizes. Applying suc h an analysis to 
the a - Centa uri system (ag ~ 18 AU), Thcbau lt et al.1 
(|2008l I2009D showed that "accretion-friendly" zones ex- 
ist only for regions within ~ 0.5 — 0.75 AU of an in- 
dividual star within the binary system. Furthermore, 
the situation would become much more unfavorable for 
planetesimal accretion if the eccentricity and precision 
of the gas disk are c onsidered (jPaardekooper et aL1l2008l : 
iBeauge et al.ll2010l) . Put simply, the planetesimals-to- 
embryos stage is strongly affected by binary effects for a 
wide range of parameter space in a circumprimary pro- 
toplanctary disk, and it is probably the main obstacle 
to the formation of planets in close binaries. Some of 
the deleterious effects of this gas drag could be miti- 
gated if the effects of a dissipa ting gas disk are taken 
into account (|Xie fc Zhoull2008D . but the effect is rather 
inefficient for small bodies and requires a long timespan 
(~ 10 6 years) for the planetesimals to reach a low Ay 
state. One possible s olution to this critical probl e m has 
been investiga t ed bv [Paardekooper fc Leinhardtl ()2010l ) 
and IXie et all (|2010b[ ). Both have shown that dust ac- 
cretion onto planetesimals may provide a safe way for 
bodies to grow through the problematic 1 to 50 km size 
range for which the perturbed environment of the binary 
can prevent mutual accretion of planetesimals. However, 
several problems still remain unclear, such as the plan- 
etesimal birth rate and the efficiency of dust accretion. 

Much of the previously mentioned work had considered 
only systems in w hich the bina r y and planetesimal disk 
were coplanar, so IXie fc Zhoul (j2009[ ) went on to con- 
sider the effects of small relative inclinations (is < 10°) 
between the binary orbit and the plane of the protoplan- 
etary disk. Overall, they found that small relative incli- 
nations signi ficantly favor accr etion through a reduction 
in AV. Then IXie et al.l (|2010af) demonstrated that if the 
disk has a small inclination, the accretion friendly zone 



can be pushed out to ~ 1 — 2 AU for a system such as 
q -Centauri. 

iMarzari et al.l ([2009D investigated cases of much higher 
inclination, a situation which may be particularly com- 
mon amongst intermediate separation binaries (40 — 100 
AU), where one does not expect there to b e sign i ficant 
alignment between disk and binary plane (lHalel Il994t 
Uensen et all 120041: iMonin et al.1 [20061 [20071) . In their 
investigation they neglected gas drag, arguing that the 
planetesimal spends the majority of it's orbit outside of 
the gas disk and hence drag effects are negligible. How- 
ever, we feel that this argument is specious, as it ignores 
two facts which are crucial. One is that the planetesi- 
mals will cross the gas disk at high inclination, i, and 
very high relative velocities, V re i, thus enhancing gas 
drag, Fd- This enhancement could be very significant 
despite the planctesimal spending less of its orbital time 
in the disk (ti n ). This is because Fd is proportional to 
the square of i 0, while ti n is inversely proportional to 
the first power of i. Thus, in total, the planetesimal 
experiences an increasing momentum change {Fd * U n ) 
with increasing i from gas drag. The other pertinent 
fact is that planetesimals in sufficiently inclined orbits 
are li kely to experience significant Kozai-driving ([Kozail 
119621 ) from the binary companion, and hence spend some 
significant portion of their time with extremely high ec- 
centricities and associated low periccnters. Given some 
typical p ower- law density profile for a proto-planetar y 
gas disk (|Havashilll98H: IPringld H98l llda fc Linll200l . 
this will mean that such planetesimals will be interacting 
with an extremely dense region of gas at pericenter, and 
hence will be likely to suffer significant drag even during 
their brief passage through the gas disk. 

Given the demonstrated importance of inclination 
and gas drag effects i n low inclination binary systems 
([Thebault et al.1 120061: IXie fc Zhoul I2009D . we feel that a 
full investigation of planetesimal evolution in highly in- 
clined binary systems requires an adequate treatment of 
gas damping. As such, we perform N-body simulations 
of planetesimal dynamics in highly inclined binary sys- 
tems, adding drag forces to simulate the effects of the 
planetesimals interacting with a gas disk. 

We structure our investigation of the planetesimal dy- 
namics in such highly inclined systems as follows: In 
Sj2] we investigate the effects of different gas disk models 
(circular and eccentric, precessing and stationary) and 
justify our choice of the simplest gas disk model. In 
we look at the behavior of individual particles as they in- 
teract with the gas disk and inclined binary companion, 
using both analytical and numerical methods. Then in 
fJH we investigate the bulk properties of a disk of plan- 
etesimals in such a highly inclined system, showing the 
overall changes in surface density that can take place, 
and their dependence on parameters such as the binary 
configuration, the radial size of the planetesimals and 
the gas disk density. Discussions of the implications for 
planet formation and the limitations of our approach are 
made in Finally, we summarize in iJB] 



' Taking the approximation for high i cases i.e., V re ; ~ iVj,, 
where Vfc denotes the Kcplcrian velocity, we get oc i 2 from 
Eqn[6] 
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2. MODEL 

FigJTJis a sketch of the orbital configuration of the bi- 
nary stars and the circumprimary disk studied in this 
paper. The companion star acts as an outer perturber 
on an inclined orbit with mutual inclination of %b around 
the inner circumprimary disk. Planctcsimals are initially 
embedded in the gas disk and are always treated as test 
particles, i.e. their growth and mutual gravitational in- 
teractions are not included and their motions are affected 
only by the gravitational forces from the two stars and 
the effects of the gas disk. The gas disk may have two 
effects on the motion of planetcsimals. One is the gravity 
of the gas disk, which is ignored in this paper and will be 
addressed in detail in our forthcoming paper (see the dis- 
cussion in § 15. 2j) . The other effect is the gas drag, which 
is generally applied to planetesimals whenever they have 
non-zero relative velocities (V re i) to the gas. Throughout 
this paper, the disk gravity is ignored and we only focus 
on the effect of gas drag to investigate the planctcsimals' 
motions in binary star systems of different orbital con- 
figurations. We use the Mercury code with the Bulirsch- 
Stoer integrator (jChamb ers 1999) for all the simulations 
except for the one in Fi gfTBI which uses the Her mite code 
as m previous papers (iXie fc Zhoull2008L l2009h in order 
to calculate the relative velocities of planetesimals. As 
gas drag is dependent on the relative motion between 
the gas and planetesimals, we must first address the gas 
motion. 

Three possible scenarios can be envisaged for the evo- 
lution of a gaseous disk under the pertur bations of an 
inclined companion (|Larwood et all 119960 . One possi- 
bility is that the gaseous disk begins to warp and it 
is disrupted by differential precession if the disk is ex- 
tremely thin. Planetesimals would evolv e in a gas-free 
enviro nment which is the case studied by iMarzari et al.l 
(|2009|) . Another possibility is that the gaseous disk re- 
mains coherent but quickly relaxes to become aligned to 
the binary orbital plane if the disk viscosity is very high. 
Then the situation becomes close to that of the coplanar 
case or slightly i nclined case, which hav e been well stud- 
ied previously (IMarzari fc Scholll 120001: iThebault et al. 
2006 1 [20081 [20091; IPaardekooper et al.ll2008l IXie fc Zhou 



20091: IXie et al.l I2010at iBeauge et al.l 120101: 
Paardekooper fc Leinhardtl I201C ) . The third possibil 



ity which is relevant for the present paper, is that the 
gaseous disk maintains a coherent structure and pro- 
cesses as a rigid body ( i.e., Ljyi S k precesses around 
L Binary as seen in FigfT]). In addition to the above 
considerations, the inner structure of the gaseous disk 
(flow line) would be significantly modified if the binary 
orbit is eccentric (iKlevi 120001 ; IPaardekooper et al.|[200l 
IKlev fc Nelsonll2008fl . 

These effects can result in very complicated gas mo- 
tions in the disk. In order to understand how gas motion 
affects the gas drag forces on planetesimals, we first de- 
scribe four gas disk models (§H3J, then give the general 
formula for the gas drag force f§ l2.2p and then finally we 
test the effects of the gas drag resulting from these four 
gas disk models f§ l2.3[) . Our goal is to obtain a reliable 
gas drag model for the planetesimals in highly inclined 
binary systems. 

2.1. Gas Disk 



We consider four models for the gas disk. 

1. The Circular Gas Disk model (CGD), where the 
gas is in a circula r orbit in the circumpri mary mid- 
plane. Following iTakeuchi fc Li n (2002), in cylin- 
drical coordinates (r, z), we adopt a gas-density 
distribution of the form 



and a gas velocity 



V g (r,z) = V^mid 



1 K B 
1 + - — 

2 V r 



qz" 



(1) 



(2) 



where h g (r) — h (r/A\J)( q+3 ^ 2 is the scale height 
of gas disk, T4, m id is the Keplerian velocity in the 
midplane, and fd is a dimensionless scaling fac- 
tor. The subscript "0" indicates the value at 1 
AU from the disk center— the primary star. In this 
paper, we take th e Minimum Ma ss Solar Nebula 
(hereafter MMSN; lHavashil (|1981l) l as the nominal 
gas disk, where fd = 1, p = —2.75, q = —0.5, 
p g o = 1.4 x 10" 9 g.cm- 3 , and h = 4.7 x 10~ 2 AU. 
We note that the CGD model is a rather crude sim- 
plification because it completely ignores the effects 
of the companion star on the gas disk evolution. 
However, to a first approximation, the CGD is rea- 
sonable and it is also the disk model most-often 



adopted by previous studies 
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2007 
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2008 


2009). 



2. The Elliptical Gas Disk model (EGD) is the same 
as the CGD model (same density profile and same 
flow speed), except that the gas is on an ellipti- 
cal orbit, with an eccentricity of e g and precession 
rate of uj g . Rigorous values for e g and ur g require 
detailed hydro-dynamical modeling on a case-by- 
case basis, which is beyond the scope of the in- 
vestigation we wish to perform. Here we follow 
an alternative strategy and derive simplified ana- 
lytic expressions for e g and uj g , based on the avail- 
able hydro-dynamical re s ults fr om previous stud- 
ies. rPaardekoope r~et al.l (|2008T ) simulated e g and 
uj g for two sets of binaries. Of particular interest 
are the simulations they performed of a 7 Ccphci- 
like binary system with mass ratio qs = 0.234, 
orbital semi-major as = 20 AU, and eccentricity 
eg = 0.3. They found that their results for e g 
and u! g depended on the flux limitcr that they used 
in their code. For the "minmod" limiter, the disk 
was "quiet" with low eccentricity e g < 0.02 and 
no precession (uj g = 0). For the "superbee" lim- 
itcr, the disk is "excited" with e g up to 0.2 and 
obvious precession {u g 7^ 0). Here we adopt the 
"superbee" results for our EGD model (the "min- 
mod" is close to CGD). U sing Fig. 7 and Fig. 8 of 
IPaardekooper et al.l (|2008f ). we estimate 

f if r< 0.5 AU 

e g ~ <^ 0.12(r-0.5) if 0.5 < r < 2.5 AU (3) 

{ 0.24 - 0.02(r- 2.5) if r > 2.5 AU, 
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and 



Lo g ~ 27r/1500yr. 



(4) 



(personal communication with Paardckooper re- 
veals that the original data underlying these simu- 
lations has been lost due to disk death, hence the 
need for us to create an approximate reconstruction 
directly from the published paper). 

We note that, in fact, the disk will never reach a 
steady state (see the Fig. 8 of Paardekooper et al. 
2008), and thus the above estimation of e g and oj g 
should be considered as an averaged approxima- 
tion. Note also that Eqn(3] and Eqn|3] arc valid 
only for the 7 Cephci-like case, a more generic ex- 
pression for e g and uj g which includes the depen- 
dency on the binary configuration would require 
more simulation results based on hydro-dynamical 
modeling. 

The Circular Gas Disk + Nodal Precession model 
(CGD+N), has the same inner structure as the 
CGD model, but adds in the precession of the as- 
cending node of the disk plane relative to the bi- 
nary orbital plane. Nodal precession occurs when 
the disk plane has a considerable t ilted angle (is) 
relativ e to the binary orbital plane. lLarwood et al.l 
(|1996f ) found that, if the Mach number of the disk 
is not too large {Ma < 30 roughly), then the disk 
can maintains its structure and undergoes a near 
rigid body precession at a rate given by 



fir. 



15M B R d 
32M A a% 



(5) 



where Ma and Mg are the masses of the primary 
and secondary stars, fi^ is the local Keplerian rota- 
tion rate and Rd is the radial size of the truncated 
disk. For the 7 Cephei like case, we set Rd = 4 
AU according to [Artymowicz fc Lubowl (|1994[ ) and 
IPichardo et all (120051 ). 

4. The Eccentric Gas Disk + Nodal Precession model 
(EGD+N), has the same inner structure as the 
EGD model and the same nodal precession as the 
CGD+N model. 

From physical considerations, EGD+N is the most re- 
alistic disk model of the four. However, as we will see in 
§ 12.31 the four disk models all result in similar gas drag 
effects in highly inclined cases. 

2.2. Gas Drag 

A spherical particle of radius R p , moving through a gas 
with a relative velocity V re i, experiences an aerodynamic 
drag force Fd given by (e.g. lArm itagc 2010h 



C 



F d = — ^-Trp g R p V re iV reh 



(0) 



where p g is the gas density and Co is the drag coeffi- 
cient. The form of the drag coefficient depends upon the 
size of the particle compared to the mean free path A of 
molecules in the gas (A ~ 0. 27 \r/AU)" 11/ ' 4 cm for the 
MMSN). If R p < 9/4A, then 



Co = r I - 

3 \ 7T 



1/2 



Vre 



where Cs ~ Qkhg is the local sound speed. This is 
called the Epstein regime of drag. For larger particles 
the Stokes drag law is valid, and the drag coefficient can 
be expressed as a piecewise function 



r 24i?e 
C D = I 24i?e 
[0.44 



-1 

-0.6 



Re < 1 

1 < Re < 800 

Re > 800. 



(8) 



(7) 



Here Re = 2R p p g V re i/ fi is the Reynolds number, and 
fi = (8/7r) 1 / 2 p g C s A/3 is the kinematic viscosity of the 
gas. 

2.3. Model Test and Comparison 

All the tests and co mparisons are conduc t ed in the 7 
Cephei-like case as in fpaardekoope r et al.l (|2008f ) with 
mass ratio qs = 0.234, orbital semimajor as = 20 AU, 
and eccentricity eu — 0.3. Although such a binary con- 
figuration is too tight to be relevant for the highl y in- 
clined case (usually as > 30-40 AU according to iHald 
((1991 1. we still choose it as our test case because hydro- 
dynamic results are currently availa ble only for such close 
binaries ((Paardcko oper et al.ll2008l ). In the following, we 
consider two sets of tests or comparisons. 

• We first test the EGD model through comparison 
(see Fig Q?l) to the hydro-d ynamical simulations of 
iPaardekooper et alj (|2008f ). In our simulation with 
the EGD model, 100 particles with R p = 1 and 
R p = 5 km respectively were initially set in cir- 
cular orbits in the coplanar plane with semi-major 
axes from 0.5 to 5 AU. As can be seen in Figl^J the 
EGD model can roughly recover the basis features 
of Paardekooper's results, such as (1) the shape of 
the eccentricity profile of planctcsimals, and (2) the 
trend that smaller particles (R p = 1 km) behavior 
more closely matches that of the gas in the inner re- 
gions of the disk (because gas drag is more efficient 
for smaller size particles and denser gas density). 
Therefore, we believe that the EGD model, though 
rather simplified, provides a reasonable means to 
include the basic effects of gas-disk eccentricity in- 
duced by the companion's perturbations. 

• We then compare the results of the four gas disk 
models (see Fig|3]). We set three planetesimals with 
radial size R p = 5 km, and with circular orbits 
starting at 3 AU (i.e., the initial semimajor axis 
a P Q~3 AU) in the gas disk plane. We consider a 
near-coplanar case (ib = 10°), an intermediately 
inclined case (is = 30°) and a highly inclined case 
(ib = 60°) for which the Kozai effect should be 
apparent. As can be scon in Figj3( only for the 
near-coplanar case does the gas eccentricity have 
any significant effect on the evolution of the plan- 
etesimal semi-major damping. That is because the 
gas drag force depends mainly on the planetesi- 
mal's orbital inclination with respect to the gas 
disk plane in significantly inclined binary cases (see 
also in Fig[5] and Fig|6]) . In addition, we see that 
the nodal precession of the gas disk has little effect 
for all the cases. This is well understood as the 
planetcsimal node always undergoes nodal preces- 
sion independent of whether the nodal precession 
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of the gas disk occurs or not. A detailed investi- 
gation of the planctcsimal dynamics, including the 
reason for large inward jumps evident in Fig|3] is 
provided in §[3] 

In summary, through the above comparisons, we see 
that the most commonly adopted model, i.e., the CGD 
model, is only a rather crude simplification in the near 
coplanar cases (is < 10°), but it is a good approxima- 
tion in highly inclined cases {%b > 30). This conclu- 
sion, although obtained from a test on a specific binary 
with a 7 Cephei-like configuration, should apply to bina- 
ries with any separation since there is not any significant 
separation dependence as seen from the above test and 
comparison. The EGD+N model, though more realistic, 
is currently available only for close binaries with sepa- 
ration of 10-20 AU, while the CGD model can readily 
be applied to any case. Therefore, given our interests in 
this paper in modeling highly inclined binaries with sep- 
aration of >30-40 AU, hereafter we adopt the CGD as 
our standard gas-disk model, and focus on it as we try to 
understand the planetesimal behavior in binary systems 
with different configurations (as, es, *s) and mass ratios 
(to)- 

3. DYNAMICS OF INDIVIDUAL PL ANETESIMALS : 
INWARD JUMPING 

3.1. Different Regimes of Semimajor Axis Damping 

Before presenting detailed numerical results, we 
present in this chapter analytical considerations in order 
to understand the different regimes in which gas drag can 
operate on planetesimals in a circumprimary disk. 

The main effect of gas friction is to force an i nward drift 
of planetesimals. According to lAdachi et al.l (|1976ft . the 
average migration rate (or the semi-major axis damping 
rate) follows the approximate relation 

x + (9) 

where a, e, i are the orbital semimajor axis, eccentricity 
and inclination (relative to the disk plane rather than the 
binary orbit plane) of the planetesimal, the hat symbol 
(e.g. e and i) denotes the average value of that variable 
taken over at least one orbital period, 77 denotes the devi- 
ation of the gas rotation from the Keplcrian velocity, and 
t is a characteristic timescale depending on the gas den- 
sity and planctcsimal radial size. For a gas disk scaled 
in proportion to the MMSN, we have r\ ~ 2 x 10 -3 and 
To can be written as 

where fd is the scaling number of the disk density com- 
pared to MMSN, and the factor (i/ho) accounts for the 
decrease in the gas density in the vertical direction. 

Eq. IHl indicates that the migration rate has contri- 
butions from three sources. One is ry, which accounts 
for the fact that the pressure supported gas disk departs 
from purely Keplcrian dynamics and only depends on the 



properties of the gas disk. The other two are the plan- 
etesimal orbital eccentricity (e) and inclination (i), which 
can be excited via various mechanisms by the companion 
star, and measure the angle between planetesimal orbits 
and gas streamlines. If, for example, i 3> e and 77, then 
Eq[9] gives da/dt oc i 2 ; planetesimals undergo a much 
faster inward migration with increasing i. This first or- 
der estimate indicates that highly inclined planetesimals 
never evolve in a manner close to the gas-free case, de- 
spite the fact that they spend more than 80 — 90% of 
their orbital time in an essentially gas-free environment 
(jMarzari et al.ll2009f ). 

Fig|4] is a schematic diagram showing which source or 
mechanism dominates the inward migration. Three main 
regions, corresponding to different regimes for gas drag, 
can basically be distinguished: (1) region S, where plan- 
etesimal migration is close to that in an unperturbed, 
single star case, (2) region B where planetesimal migra- 
tion is dominated by the binary perturbations, and (3) 
region SB, a transition region between B and S. In addi- 
tion, region B can be divided into three subrcgions corre- 
sponding to three different mechanisms which can excite 
the e and/or i. The precise definitions of all these regimes 
is as follows: 

• Region S: where |e 2 < rj 2 and \i 2 < n 2 . For this 
low dynamical excitation case, according to Eq[9j 
the migration rate is dominated by r]. For MMSN- 
like gas disks, Region S corresponds to planetesi- 
mals that are very close to the primary as compared 
to the distance of the companion, i.e. a/ as <C 1. 
In this case, e and i remain very small and planetes- 
imal orbits are close to what they would be in an 
unperturbed single-star case. Note here that the 

inward migration rate (Eqn(5]) da/dt oc r/ 2 a/TQ oc 
a -1 , meaning that the inward migration is acceler- 
ated in region S. 

• Region B: e and/or i are significantly pumped up 
by the companion's perturbations and thus domi- 
nate the damping rate of planetesimal semimajor 
axis. Depending on the specific mechanism that 
pumps up e and/or i, region B is divided into three 
subregions: 

— subregion Bl: where e and/or i are excited 
by secular perturbations, but Kozai effects arc 
absent, i.e., is < 39.2°. If we assume the aver- 
age values of e and i approximate their forced 
values e/ (e/ <^ 0.1, see the Appendix) and 
if (if ~ i B < 39.2°), then region Bl corre- 
sponds to the middle part of FigH] 

— subregion B2, the Kozai regime, i.e., i ~ is > 
39.2°. In such case, both e and i can be highly 
excited via the Kozai effect. 

— subregion B3, where e is excited and it is in- 
dependent of i. This region may correspond 

to a mean motion resona nce scenar io. An 

example can be found in llda fc Lid (fl~996h . 
which shows the combination of gas drag and 
Jupiter's mean motion resonances can affect 
the structure of the asteroid belt. 
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• Intermediate region SB: where |e > 77 and > 
?y 2 but y|e 2 < 7yand|i 2 < 77. Region SB corre- 
sponds to the transition case, where e and i are only 
mildly excited, meaning that the migration rate is 
determined by the combination of 77, e and/or i. 

The dashed line in FigJH i.e., y|e 2 = |i 2 , divides the 
entire plot into two parts: In the upper-left e dominates 
the damping of semi-major axis, while in the bottom- 
right of the plot i contributes more than e. We define the 
boundary between regions Bl and B3 as being e = 0.1, 
since the maximum forced eccentricity (e/) for a typical 
binary system is ~ 0.1 (see the Appendix for detail). 

Note that Figgis based on EqJS] which is derived from 
the most simplified gas drag model, i.e. the circular gas 
disk model (CGD). According to the test in §2.3 (see 
also in Fig|3|), the CGD model is a good approximation 
only for relative large i. Therefore, hereafter we focus 
only on these large i (or i B ) cases, and simply classify 
them into two categories: Kozai-off (which corresponds 
to the right-hand section of the Bl region) and Kozai-on 
(the whole of the B2 region). We don't explore region 
B3 in the present paper, as it is likely only to be relevant 
to some isolated (i.e., in resonance) or unstable regions 
(i.e., close to the companion star) of the protoplanetary 
disk. 



3.2. Kozai-Off Regime 

In the following sections, we turn to numerical simula- 
tions using the circular gas disk model (CGD) which was 
tested and justified in §[3] 

We first define our standard system parameters upon 
which we base most of our variations. We take Ma = 
2Mb = A/©, a B = 50 AU, e B = 0.3, i B = 30°, 
fd = 1.0, and R p = 5 km. The mass ratio and eccen- 
tricity are chosen as the typical values for known bi- 
naries ([Duauennov fc Mayor! Il99ll ) . The separation is 
chosen as it is wide enough to be relevant for large incli- 
nation and tight enough to have the binary companion 
cause considerable perturbations in the 0.1 — 10 AU re- 
gion, where the core-accretion model of planet formaiton 
prefers. i B = 30° is chosen as it is large enough for the 
validity of CGD model as tested in § 12.31 but not so large 
as to induce the Kozai effect. Such a binary configura- 
tion leads to a stable orbit limit of ~ 1 AU around the 
primary star dHolman fc Wiegert]|1999D . 

This high-i case corresponds to the right-hand side of 
the Bl region in FigUJ FigJS] presents the result of a 
numerical integration showing the typical behavior of a 
planetesimal orbit in this regime (see also in the middle 
panel of Fig(3} but note different binary parameters). As 
expected from linear secular perturbation theory, e and 
i undergo oscillations with a period equal to the secu- 
lar perturbation time scale, i.e. t sp . The value of t sp 
derived from Fig(5] is in good agreement with the fir st 
order approximation derived bv lThebault et all (|2006D : 



1/2 /MbV 1 

yr, 



Ma 



Mb 

M q 



(11) 



where Ma, M b , and Mq are the masses of the primary, 
secondary and the Sun respectively. As the relative ve- 
locity V re i between the planetesimal and g as can be ap- 
proxim ately estimated as (see Eqn 4.20 of lAdachi et al.1 

MM) 
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this means that in the regime where i 3> e(and?7), V re i 
and the gas drag force Fd (EqnJB]) will also follow the 
same periodic oscillation as i. This oscillation of gas 
drag damping causes the planetesimal semimajor axis to 
damp stcp-by-step as seen in the bottom panel of FigJS] 
In addition, as can be seen in FigJS] the oscillation 
amplitudes of i an d e are under-dampin g because of gas 
drag. According to lAdachi et al.l (11976D . the damping of 
i and e follows, 



TH— r^2 — — 

e dt i dt 



5-2 . 1"2 
8 6 + 2 J 



if 



1/2 



(13) 



In the present case, where i>e (and 7/), the timcscale of 
i damping can thus be approximated by 



^ dp 



p3/2^£ 



(14) 



di/dt i 
Substituting Eq|10]into EqlH t dp can then be rewritten 



-1/2 



yr (15) 



Depending on whether the secular perturbation or the 
gas drag damping dominates, FigJS] can be divided into 
two parts with the boundary between these two regions 
being roughly set by the last turn-over point of the mi- 
gration curves: E.g. for the red-curve in FigjSJ the coor- 
dinates of the critical time and semi-major axis at which 
this occurs is t c ~ 5.75 x 10 5 yr and a c ~ 1.5 AU for 
a p0 = 7 AU. For times earlier than t c or semi-major 

axes beyond a c , i is periodically excited, and thus the 
planetesimal periodically undergoes a step-wise inward 
migration with the steep drops in semi-major axis cor- 
responding to the peaks in i. For times greater than 
t c or semi-major axes inside a c , i is smoothly damped 
to a low value, thus the planetesimal migrates inward 
smoothly and slowly. 

We can also estimate a c by equating t sp and td p fEqfTTI 
and Eq fTB]) . which gives a "zeroth-order" estimation of 
a c o as, 



a c o ~ 0.9 x f 2 d /9 ( 



km 



as 

10AU 

-2/9 



2/3 



(l-e|)V3 



Vf (M± 



\M E 



2/9 



AU. 



(16) 



We call a c o "zeroth-order" because this timescale analy- 
sis gives the correct dependency of a c on the various pa- 
rameters (for example, a c is independent of a p o in FigJS]), 
but not the correct absolute coefficient. Through many 
sets of simulations, such as Fig[5]but with different stellar 
mass ratios (Ma/M b =0.2, 0.5, 1.0), binary semimajor 
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axies (a B =10, 20, 50, 100, 200 AU), binary eccentricities 
(e B =0, 0.3, 0.6), binary inclinations (i B = 20°, 30°, 39° 
and 150°), gas disk densites (fd =0.05, 0.5, 1), and plan- 
etesimal radial sizes (R p =1, 5, 100 km), we found similar 
results for both the prograde and retrograde cases (see 
the bottom panel of Figj5|) and the empirically measured 
form of a c is just systematically smaller by ~ 33% than 
a c0 , i.e., 

2 

(ir 8 (^r »■ ™ 

We note that t dp is larger by just about one order of 
magnitude than t sp at a c . 

The physical meaning of a c is that it indicates the turn- 
over point where the planetesimal semi-major axis tran- 
sits from the excited (usually fast) damping regime to the 
quiet (usually slow) damping regime (also corresponding 
to the boundary between the region Bl and SB). The 
importance of a c lies in the fact that it determines the 
location in the proto-planetary disk where planctcsimals 
will pile up (see §2] for details). 

3.3. Kozai-On regime 

The standard setup for this case (corresponding to the 
B2 region in Fig|4]) is chosen as similar to that in § 13. 2i 
but with i B = 60° instead. As i B is larger than the 
critical value 39.2°V(R), the Kozai effect is turned on 
(Kozai 1962), and e can be excited to a degree which is 
comparable with i. In such a case, both e and i have 
major contributions to the planetesimal migration rate 
(see Eqn[9]). 

The most interesting result is that large eccentricity 
oscillations can lead to very fast inward migration of the 
planetesimals towards the inner disk. As can be seen 
in Fig|51 a planetesimal starting at a p o = 3 AU spirals 
down to 1 AU in ~ 10 4 years. For the a p o = 5 AU case, 
the planetesimal starts closer to the companion star and 
is subjected to much stronger Kozai perturbations. Ec- 
centricities can quickly reach values as high as 0.7 that 
places the planetesimal periccnter in inner gas rich re- 
gions where it will in addition reach very high V re i values. 
This leads to a very fast, almost "jump" -wise inward mi- 
gration from 5 AU to - 0.35 AU. For the a p0 = 8 AU 
case, the eccentricity also reaches 0.7 but, due to the 
larger initial semi-major axis, the planetesimal cannot 
directly reach very dense inner regions and thus first un- 
dergoes several small jumps, until finally making a bigger 
jump to ~ 1 AU. 

All the inward jumps occur when e and/or i is excited, 
leading to a surge in the experienced gas drag force (F d ). 
Although the e contribution to Fd is significant, it seems 
that i is still the major component of Fd. Note i is with 
respect to the disk plane rather than to the binary orbital 
plane, and thus the periodic variation of i shown in Fig|5] 
is caused by the Nodal precession of the planetesimal 
orbit with a period of two times that of the Kozai cycle 
period. 



For all cases, once the planetesimals complete their 
jump, they will approximately stabilize at that semi- 
major axis, with very small residual e and i and thus 
only very slow inward migration. They no longer suffer 
from the Kozai effects because they have jumped into 
an inner region where the gas is so dense that gas drag 
damping dominates over Kozai excitation. Unlike the 
Kozai-off case, the turn-over semi-major axis, where the 
planetesimal stops jumping, is no longer independent of 
a p o and not equal to a c (Eqnll7|) . FigJT] shows how it 
varies with a p o in the Kozai-on case. 

• If a P Q < 1.5 AU, the turnover semi- major axis 
is around a p0 , which means that no "inward- 
jumping" occurs to the planetesimals in the very 
inner region; Kozai effects is fully turned off be- 
cause such dense gas and hence strong gas drag 
damping dominates over Kozai excitation. 

• If 1.5 < a p o < 8 AU, the turn-over semi-major axis 
follows a characteristic oscillatory pattern, whose 
amplitude varies between a lower limit at a c a ow \ ~ 
0.3 AU and an upper limit at a c ~ 1.5 AU roughly; 

• If Q-pO It 8 AU, the planetesimal is initially located 
close to the boundary of the stable orbit (~ 10 AU), 
thus the results become chaotic. However most, 
if not all, of the turn-over points remain confined 
within the range between a c (i ow ) and a c , i.e., 0.3- 
1.5 AU. 

In FiglSl we study the dependence of the lower limit 
of the turn-over semi-major axis (cl c (Iow)) on a B and i B . 
Empirically we found that 

f0.50a c if i B = 50° 
a c n ow) ~ <^ 0.20a c if i B = 60° (18) 
[0.10o c if i B = 70° 

where a c is given by Eqn ll7l For even larger i B {7Q° — 
90°), the a c n ow ) turns to be larger. The a c (j ow ) for 
prograde and retrograde cases are symmetrical about 
i B = 90° . The combination of a c (; ou ,) and a c determines 
a region for the Kozai-On case, into which the planetes- 
imals are likely to jump and pile up. 

4. BEHAVIOR OF A SWARM OF 
PLANETESIMALS: PILING UP 

4.1. Initial Set-up 

In §[3j we focused solely on the dynamics of individual 
planctcsimals. In this section, we investigate the behav- 
ior of a swarm of planetesimals, focusing on the evolution 
of the planetesimal disk's surface density S, under the 
coupled effects of gas drag and binary perturbations. 

In each simulation, all planctcsimals are initially set 
onto a near-coplanar circlular orbit around the primary 
star with orbital eccentricity and inclination (relative to 
the gas disk plane) < e = 2i < 10~ 4 . All other angular 
elements are randomly distributed. The planetesimals' 
scmimajor axes are uniformly distributed in a range ,i.e., 
din < a < a out with the same radial size (R p ). We set 
a-in = 0.05 AU, within which planetesimals are removed 
from the simulation. a ou t is set as the boundary o f the 
stable orbit according to iHolman fc WiegertJ (|1999f ). In 
order to obtain reliable statistical results, we initially put 
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Npo = 100 planctcsimals in every semimajor axis bin of 1 
AU width. By monitoring the planctcsimal number (N p ) 
in each bin, we can estimate the change of planctcsimal 
surface density, i.e., E/Eo- 

Note that during the statistical process (i.e., count- 
ing N p ), we asign two weighting factors (wi, IU2) to 
each planetesimal (the tracked test particles), where 

— 1/2 

W\ = cipo leading to an initial surface density profile 

of S oc a~ 3/2 as in MMSN, and w 2 = 4.2 if a pQ > 2.7 
AU otherwise w% = 1, accounting for the solid density 
enhancement beyond the snow line. The gas disk model 
is the CGD model as described in § 12.11 

We conduct a large set of simulations in which 
we vary the stellar mass ratios (q = M B /Ma = 
0.2,0.5,1.0), binary orbital semi-major axes 
(a B = 10,20,50,100,200 AU), binary orbital eccentrici- 
ties (e B = 0,0.3,0.6), binary orbital inclinations [i B = 
10°, 30°, 50°, 60°, 70°, 80°, 90°, 100°, 110°, 120°, 130°, 150°), 
and the planetesimal radial sizes (R p = 5,100 km). 
For the convenience of describing the results, we first 
define a standard case, then we vary one of parameters 
of the standard case each time to see the effects of the 
parameter on the result. 

4.2. the Standard Case 

For the parameters of the standard case, wc choose q = 
0.5, a B = 50 AU, e B = 0.3, i B = 60°, and R p = 5 km. 
The results of the standard case (highlighted with thick 
blue axes) are shown in the second panel of the right- 
hand column of Figj9]and in the middle panels of FigflOl 
and FigfTTJ the results of which can be summarized as 
follows. 

• In the first 2 x 10 4 yr, the planetesimal disk profile 
varies a lot. In the outer region the profile be- 
comes very noisy with modest fluctuation around 
So- This complicated structure is due to the step- 
wise jumping as shown in Fig|5] (see the bottom- 
right panel) and the chaotic behavior near the outer 
boundary (see Fig(7J) . The modest density surge 
near ~ 2.7 AU is due to the initial density jump at 
the "snowline" ai ce = 2.7 AU. In contrast, in the 
inner region (~ 0.3 AU) a significant density en- 
hancement (up to E/Eo ~ 10) is observed, caused 
by the planetesimal jumping shown in Fig|5] (see 
the bottom-middle panel) and Fig(7] 

• At times ~ 10 5 yr, the planetesimal disk appears 
truncated near a c (EqfT7|. Within ~ a c , the 
entire surface density is enhanced by 1 order of 
magnitude, i.e., E/Eo ~ 10, while beyond ~ a c , 
the surface density is entirely depleted down to 
E/E < 0.1. This result is expected from Fig 13 
which shows that most, if not all, planctcsimals 
jump into the inner region within a c . 

• At later times (t > 10 5 yr), the planetesimal disk 
seems to reach a quasi-cquilibrium state: the sur- 
face density maintains a (daleth) shape profile 
with the sharp outer boundary from ~ a c slowly 
moving inward. This "orderly" evolution is ex- 
pected because, as show in Figj6] ( §[3T3j , once the 
planctcsimal has completed a jump within a c , its 



e and/or i would no longer be excited by the bi- 
nary perturbations and thus it undergoes slow and 
smooth inward migration. In other words, the 
planetesimal's inward migration has transited from 
Bl regime to SB and finally to the S regime as 
shown in Figf?] (§ I3.1[) . As the inward migration 
rate is inverse proportional to its semimajor axis 
(see § 13. ip in the S regime, the surface density of 
planetesimals will then decrease from inside to out- 
side on a long timescalc (depending on the gas den- 
sity and the planetesimal size). 

4.3. Parameter Exploration 

Through comparisons to the standard case, we will dis- 
cuss the effects of the modeling parameters, such as i B , 
a B , e B , q B , R p and f d - 

4.3.1. Effects of i B 

We vary the parameter i B of the standard case from 
10° to 150° but hold all other parameters constant. Our 
results are plotted in Fig[9] (the standard case is in the 
second panel of the right-hand column of Fig{9j outlined 
in bold), and can be summarized as follows. 

• Kozai-On case (40° < i B < 140°). Similar to 
the standard case, the planetesimal disk reaches 
a quasi-equilibrium state with a (daleth) pro- 
file on a timescale of ~ 10 5 yr (This process is 
a little bit quicker for larger i B cases). In all 
cases, the planetesimals pile up within a c , leading 
to an average surface density enhancement up to 
one order of magnitude, i.e., E/Eo ~ 10. Note, re- 
sults of the prograde cases (i B = 50°, 60°, 70°, 80°) 
and their corresponding retrograde cases (i B = 
130°, 120°, 110°, 120°) are very similar (see also 
in FigJH]). The only slight difference lies in the 
marginally quicker clearing of the outer disk. 

• Kozai-Off case (i B = 30°, 150°). In contrast to 
the Kozai-on case, the planetesimal disk seems to 
evolve to a unimodal profile on a timescale of ~ 
10 5 yr. This unimodal profile is maintained with 
a peak enhancement (E/Eo ~ 10) at ~ a c , slowly 
moving inward over time. Another feature which 
distinguishes this from the Kozai-On case is that 
the disk will not be significantly truncated with 
respect to the initial disk. Note again, the prograde 
case (i B = 30°) and its corresponding retrograde 
case (i B = 150°) have very similar results. 

• Near-coplanar case (i B = 10°). Only modest 
density enhancement occurs in ~1.5-2.7 AU be- 
cause planetesimals pass through the "snowline" 
(a ice = 2.7 AU). 

4.3.2. Effects of a B 

We vary the parameter as of the standard case from 
10 to 200 AU but hold all other parameters constant. 
Results are plotted in Fig[TU] (the standard case is in the 
middle panel highlighted with blue bold axes) , which are 
summarized as the following. 

• Similar to the standard case, in all cases (except the 
a B = 200 AU case) the planetesimal disks reach a 
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quasi-cquilibrium, i.e. a (dalcth) shape surface 
density profile with all planctcsimals piling-up in a 
tight, dense inner region which is truncated near 
a c . In the a B = 200 AU case, the planetesimal 
disk is not truncated and there are still many plan- 
etesimals left beyond 20 AU after 10 6 yr. This is 
mainly due to the very low gas density there (thus 
weak gas drag) and due to the long binary separa- 
tion (thus weaker perturbations). In addition, the 
surface density enhancements in the a B = 10 and 
20 AU cases are significantly lower (E/Eo ~ 5) 
as compared to the standard case (a B = 50 AU, 
S/So ~ 10) because the initial planetesimal disks 
are much smaller (thus less total planctcsimals) in 
such close binary systems. 

• The time for the planetesimal disk to reach a quasi- 
cquilibrium (""1-likc shape) surface density profile 
(we define it as t p u e ) is very sensitive to a B . Em- 
pirically and to a first order of estimation, t p u e is 
roughly equal to the secular perturbation timescale 
at a c . Substituting the semimajor axis in Eqn llll 
with Eqn ll7i t p u e can be expressed as 



t 



pile 



x a-4) 



VAU/ 

V6 / n r \ -2/3 



(Ma\ l (Mb 



yr, 



(19) 



which depends mainly on as and weakly on e B , 
M A and M B . For a B > 200 AU, the pile timescale 
tpiie > 1 Myr, which is comparable or even longer 
than the lifetime of the gas disk. Therefore, the 
jump-piling effects would be very limited for binary 
systems with a B > 200 AU. 

4.3.3. Effects ofeB, qB, R P and fd 

The results are shown in FigfTTl with the standard case 
located in the middle panel with blue bold axes. We can 
obtain the effects of e B and q B through the comparison 
of the three middle vertical panels and the comparison 
of the three horizontal panels, respectively. In addition, 
we can compare the cases of R p = 100 km and R p = 1 
km (the two bottom-corner panels) to the standard case 
(the middle panel with blue bold axes) for the effects of 
R p . Since increasing the planetesimal radial size (R p ) 
is equivalent to decreasing the gas density^) from a 
dynamical view, thus we obtain the effect of fd at the 
same time. 

• Effects of e B (q B ). Increasing e B (q B ) leads to 
quicker evolution of the disk density profile (shorter 
tpiie as also seen the dependence in Eqn riU| 
and slightly lower surface density enhancement 
(E/Eo) due to the stronger binary perturbation 
and smaller initial planetesimal disk size. 

• Effects of R p (fd)- Increasing R p (decreasing fd) is 
a double-edged sword when it comes to the plan- 
etesimal pile-up. On the one hand, it directly re- 
duces the gas drag acceleration on the planctcsi- 
mals and thus slows their "jumping processes" and 
inhibits the planetesimal piling (see the leftover 



case in Fig|TT]). On the other hand, once planetesi- 
mals have finished their jump into the inner region, 
they will experience much stronger braking of their 
further inward migration because of the larger R p 
(smaller fd), which favors the planetesimal pile- up 
(see the nearly fixed boundary around a c , and the 
progressive increase of the surface density within 
a c in the R p = 100 km case). The opposite trend 
can be found by decreasing R p or increasing fd (see 
the case of R p = 1 km) . 

4.4. Analytic Estimate o/E/Eo 

As we have shown above, the inward migration or 
jumping process always lead to the formation of a com- 
pact, dense disk with an outer boundary around a c . 
Therefore, the surface density enhancement in the in- 
ner region can be estimated by assuming that all of the 
material initially having a > a c is transported to a < a c , 
leading to an enhancement factor 



E f a Qc 27rrE dr' 



(20) 



where ai n (here set as 0.1 AU) is the inner boundary of 
the planetesi mal disk and ad is t h e init ial outer bound- 
ary following iHolman k Wiegertl (|1999l ). Considering a 
power law disk profile as Eo oc r Q (o / 2), then Eql20l 
can be rewritten as 



E 
Eo 



if a c < a d < ai 



if a c < a ice < a d , (21) 



if a ice < a c < a d , 



material in the outer disk of the R v 



100 km 



where ai ce is the location of snow line and /j ce is the solid 
density enhancem ent beyond the snow line. Following 
llda k Lin! (|2004ft . we take a lce = 2.7 AU and f tce = 4.2. 

Fig 112] shows the results of Eqj2TJ We see that the 
analytical estimates are consistent with the N-body sim- 
ulation results by comparing the red solid curve of Fig ll2l 
to FiglTUI Piling effects are most significant for binary 
systems with a B ~ 100 AU and for planetesimal disks 
with flatter initial disk profiles. While the profile we con- 
sider in detail (a = 3/2) have surface density enhance- 
ments ~ 10, for disk profiles flatter than Eo oc r -3 / 4 , the 
surface density enhancement can be around 2 orders of 
magnitude. 

5. DISCUSSIONS 

5.1. Implications for Planet Formation in Highly 
Inclined Systems 

In this paper, we have shown, for the highly inclined 
(*s 30°) binary systems, that planetesimals from the 
outer disk can "jump" inwards and pile up in the inner 
disk. The typical result of this jumping-piling process 
is the formation of a smaller and denser planetesimal 
disk around the primary, with the truncation boundary 
near a c (Eqnll7[). within which the surface density en- 
hancement is up to E/Eo ~ 10. This significant change 
in the planetesimal disk will definitely affect the subse- 
quent planet-formation processes— planetesimal growth, 
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and the final composition and configuration of the plan- 
etary system. 

The two most important parameters which govern 
planetesimal growth are, first, the relative velocities 
(V re i) among planetesimals, and second, their colli- 
sion rate. The former determines the outcome of 
planetesimal-planetesimal collisions, i.e., accretion or 
erosion; while the latter controls the speed of plan- 
etesimal accretion and/or erosion. We investigate in 
Fig[l3] the evolution of these two parameters during the 
jumping-piling process for the standard binary config- 
uration (Ma = Mq, qB = Mb/Ma = 0.5, a B = 50 
AU, e B = 0.3, i B = 60°). In Figdl 10,000 planetesi- 
mals with a centered Gaussian size distribution between 
R p = 1 km and R p = 10 km, are initially randomly dis- 
tributed between 0.3 AU and 10 AU (flat distribution) 
with near circular and coplanar orbits. As can be seen, 
in the first ~ 10 4 yr, the average V re i has a steep increase 
up to 10 3 ms^ 1 because of the strong binary perturba- 
tions, and the collision rate is relatively low because of 
the large spread of planetesimals at the beginning. After 
this point (t ~ 10 4 yr), planetesimals begin to jump in- 
wards and pile up in the inner region (within ~ 1.5 AU), 
leading to a surface density enhancement of up to one 
order of magnitude, a collision rate enhancement of up 
to three orders of magnitude, and a reduction of aver- 
age V re i down to 40 — 50 ms _1 . If we adopt a velocity 
of 50 — 100 ms" 1 as the thresho ld for planetesimal ero- 
sion as in (jMarzari et al.l I2009D , then we see that the 
jumping-piling process moves most planetesimals from 
an erosional regime to an accretional regime. For this 
reason, we expect fast planetesimal growth to occur in 
the inner region. 

In fact, everything (such as V re i, e, and i of planetes- 
imals) within a c , is getting close to that in single star 
cases. For example, after planetesimals have completed 
the jumping process as shown in Figj6l their eccentric- 
ities (e) are all below 10 -3 , which are comparable to 
that in single star cases (see Eqn.(10) of IKokubo fcTdal 
(|2000D ) and their inward drifts transition to the SB or 
S regime. Therefore, the process of planetesimal growth 
within a c would probably follow a similar way as in sin- 
gle star systems. Such a process leads eventually to the 
formation of planetary embryos (via the runaway and oli- 
garchic growth processes), resulting in a number of em- 
bryos which have isolated themselves from one another 
due to the accretion of all the planetesimals in an annular 
feeding zone aro und themselves. This embryo isolation 
mass is given by (llda k Linll200l 



A/ iso ~ O.W/i 



3/2 f 3/2 



So 



3/2 



3/4 



M®. (22) 



where fi ce ~ 4.2 denotes the solid density enhancement 
beyond the ice line a,i ce = 2.7 AU, and E/E denotes the 
solid surface density enhancement relative to MMSN. 

Clearly, in a MMSN, at 1 AU the isolation mass is 
0.16M®, leading to the idea that in the Solar System 
there may have been a chain of such ^Mars-mass em- 
bryos in the terrestrial planet region, which subsequently 
went on to collisionally evolve in a chaotic growth phase. 

However, in the binary scenario we are considering, 
we have seen typical surface density enhancements of 
S/S ~ 10 (can be larger for flatter disk profiles as 



shown in Fig fT2"l) . leading to an isolation mass at 1 AU 
~5x /d 3 / 2 M®. Moreover, the initial surface density of 
the disk could be above that of the MMSN, so it is en- 
tirely plausible to have fd = 1 — 3, giving isolation masses 
at 1 AU in the range 5 — 26 M® and typical separations 
10r H ~ 0.25 - 0.43 AU. 

From the first order estimation above, one could 
anticipate that in the highly dense inner-disk regions 
of these inclined binary systems, several giant embryos 
could form. If these cores grew before the depletion of 
the gas disk, then several gas-giants could start to form, 
otherwise, a series of stalled super carths/ncptuncs 
may form. In addition, as the planetesimal disk is 
severely truncated after the jumping-piling process, it 
is highly probable that such planets would be born 
in a dynamically compact region with a short sta- 
bility timescale (|Chatteriee et al.1 [2008). After the 
depletion of the gas disk, the final configuration of 
the planetary system would then be further shaped 
by dynamical ev olution mechanisms such as planet- 
planet scattering (iRasio k Fordl lBM iFord et all 120011: 
Malmberg et all 120071), Kozai migration (IWu k Murrav| 



2001 iFabrvckv k Tremaind 120071: iTakeda et al 
2009) and secular chaos (iMichtchenko et all 120061 



lithwick k Wull2Cltnt IWu k lithwicldl20T0ir 



5.2. Limitations 

First, as presented in §[2l we adopt a very simplified 
gas disk model (CGD) where the gas is in a circular or- 
bit in the circumprimary mid-plane without feeling the 
companion star as it is in a single star system. In real- 
ity, the gas disk structure should be significantly modi- 
fied by the binary perturbations. Nevertheless, we found 
that the planetesimal dynamical behaviors are very sim- 
ilar under different models of disk structure for highly 
inclined (30° < is < 150°) binary systems . 

Second, the growth and/or erosion of planetesimals due 
to mutual collisions are not considered in the present pa- 
per. Beyond a c , planetesimals arc excited onto orbit with 
large inclinations, leading to their collisional timescale 
becoming longer than 10 6 yr (according to Eqn. 12 of 
Xie et al. 2010), which is larger than the piling timescale 
tpiie (Eqn [T^|) . Therefore, outer planetesimals can jump 
and pile up in the inner region before they collide with 
one another. On the other hand, in the piling region 
(< a c ), planetesimals may grow up quickly if collisions 
are considered because of the very high collision rate and 
low relative velocities as shown in Fig ll3l As they grow 
up, their inward migrations slow down, which can make 
stronger piling effect in the inner region. 

Third, in order to estimate the degree of planetesi- 
mal piling effects, we have assumed that all planetesi- 
mals have already been born before the jumping-piling 
effect takes place. This assumption is reasonable only if 
the planetesimals birth rate is very high over the whole 
disk. However, some theoretica l and observational works 
(ICuzzi k Weidenschillingll200& fScottl f2007t iNatta et all 
120071 : iChauvin et all l2010h suggest that such a transi- 
tion process, from dust to planetesimal, may take as 
long as a few 10 6 years, which is comparable to or even 
larger than the piling timescale t p u e fEqn fTT))) . In such 
a case, much smaller, maybe m m-sized, particles piling- 
up through other mechanisms (jYoudin k Chiang! 120041) 
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should be considered. 

Finally, perhaps the main uncertainty in the planetes- 
imal "jumping-piling" effect comes from not including 
the effect of gas gravity on the planetesimals, i.e. we 
have neglected disk self-gravity. One possible effect of 
gas gravity might be that it can pull highly inclined plan- 
etesimals back towards the mid-plane of the gas disk, so 
that the Kozai effect (and hence the "jumping -piling" ef- 
fect) might be inhibited (|Fragner et al.l 120111) . However 
this effect usually happens for massive gas disk. For the 
standard case we are considering in this paper (ag = 50 
AU), the disk mass should be larger than ~ 6 MMSN 
to switch off the Ko zai effect as shown in the Fig. 16 of 
iFragner etaTl (|2011l ). For a nominal disk with 1 MMSN, 
we thus expect that the "jumping-piling" process will 
still be efficient even if the disk gravity is included. In 
fact, a jumping trend of p lanetesimal semimajo r axis can 
be discerned in Fig. 11 of IFragner et al.l (| 2 1 If ) . In addi- 
tion, as is effectively shown in the bottom-right panel of 
FigQTJin our paper, the "jumping-piling" process is still 
efficient even for a very tenuous gas disk of less than 0.05 
MMSN. Therefore, it is plausible that the planetesimal 
"jumping-piling" process described in this paper would 
be most important during the gas dissipation phase. 

6. SUMMARY 

In this paper, we have investigated the intermediate 
stage of planet formation - from planctcsimals to plane- 
tary embryos - in highly inclined binary systems with a 
focus on the effects of gas drag on planetesimal orbits, 
especially the evolution in the semimajor axis distribu- 
tion. 

First, we justify our choice of a simplified circular gas 
disk (CGD) for modeling gas drag force in highly in- 
clined cases. Through numerical tests and comparisons 
(see Figf2] and Figf3]), we have shown that eccentricity 
and pericenter precession of the gas disk are not impor- 
tant to the effects of gas drag on a planetesimal for bi- 
nary systems with inclination 30° < i B < 150°. The 
major factor in determining the magnitude of the gas 
drag forces is the mutual inclination i between the plan- 
etesimal and the gas disk, which is mainly controlled by 
i B and is independent of the specific disk model. 

Next, we analyze the dynamics of individual planctcsi- 
mals under the influence of gas drag which leads to the in- 
ward migration of plantcsimals. Depending on the domi- 
nant mechanism which drives the planetesimal eccentric- 
ities and inclinations, very different migration behaviors 
can take place. There are three main regimes to consider 
(see Figf?]): Regime S, where the situation is close to that 
in single star systems; Regime B where inward migration 
is dominated by the binary's perturbations, and an in- 
termediate transition regime SB. In addition, regime B 
can be divided into three sub-regimes, two of which, the 
Kozai-off regime and the Kozai-on regime, are studied in 
detail in the present paper. 

For the Kozai-off regime (i B < 39.2° or i B > 140.8°), 
planetesimal inward migration is determined only by the 
excitation and evolution of its orbital inclination (i). 



There is a turnover semimajor axis (a c , see Ean ll7l and 
Fig©: where planctcsimals transit from a fast inward mi- 
gration (regime B) to a slow one (regime SB or S). The 
value of a c is independent of planetesimal initial semima- 
jor axis (a p o), and it can be estimated by equating the 
timcscale of secular perturbation (t sp , see Eqn lllj) and 
the timescale of gas drag damping (td p , see Eqn [T5)l . 

For the Kozai-on regime (39.2° < i B < 140.8°), plan- 
etesimals can be excited to orbits with high eccentric- 
ities and high inclinations (or even become retrograde 
with respect to the gas disk). In such a case, planetesi- 
mal inward migration can be very sudden. Planetesimals 
in the outer disk can directly jump into the inner disk 
where they then slow their inward migration as their or- 
bits circularize and become coplanar due to strong gas 
drag damping (see FigEJ) . Unlike the Kozai-off case, the 
turnover semimajor axis for this migration transition is 
no longer independent of a p o, but it is always located in 
a range between a c ^ ow ^ (see, Fig[7|and FigfS]) and a c . 

We then studied the behavior of a swarm of plan- 
etesimals for both Kozai-on and Kozai-off cases with an 
exploration of many parameters (see Figf9j FigQUl and 
FigfTTj), including the binary mass ratio (q B ), semimajor 
axis (as), eccentricity (e B ), inclination (is), gas density 
scale (f c i) and planetesimal radial size (i? p ). A robust re- 
sult is that planctcsimals migrate/jump inwards and pile 
up in the inner region, leading to a smaller and denser 
planetesimal disk with a truncatation boundary near a c 
and a surface density enhancement or order X/Xo ~ 10. 
The timescale for such a "jumping-piling" process to op- 
erate, tpu e (see Ean fTU)) . is mainly dependent onos. 

Applying such "jumping-piling" effects to planet for- 
mation, we expect a growth-favored region for planetes- 
imals in the inner disk (within a c ), where collision rates 
are high and relative velocities are low (see FigfT3|). Such 
conditions may lead to the formation of embryos suffi- 
ciently massive to undergo runaway gas accretion (pro- 
vided they form prior to the dissipation of the gas disk) . 

This work, focusing only on the effect of gas drag on 
a planetesimal, is our first step towards understanding 
planet formation in highly inclined binary systems. Fu- 
ture work needs to account for important physical factors 
such as the effects of gas disk gravity, planetesimal accre- 
tion/fragmentation, and the time at which planetesimals 
emerge and the jump-piling process begins. 
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Fig. 1. — Sketch of the orbital configuration of the highly inclined binary system studied in this paper. Loi s fc and L Binary are the 
normals of the gas disk plane and the binary orbital plane, respectively. The angle is between Lr)i a k and L Binary corresponds to the 
mutual inclination of the two planes. 
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Fig. 2. — Comparison between the results given by the hydrodynamical model (top: Paardckoopcr ct al. (2008)) and the approximate 
EGD model (bottom). We display snapshots of planetesimal eccentricity vs. semimajor at t=1300 and 2300 yr. The black solid lines 
indicate the forced eccentricity ef = (5aes)/(4as ). The black dashed lines show the eccentricity distribution of the gas. We see that the 
results from the EGD model can recover the basic features of the results by Paardckoopcr ct al. (2008), thus providing us with a means to 
simply and quickly add the main features of the time-consuming hydrodynamical results into a fast N-body code. 
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Fig. 3. — Comparison between the results obtained from 4 different disk models (CGD, CGD+N, EGD. EGD+N) for three cases with 
is = 10°, 30° and 60°. We plot the evolution of a planetesimal's semi-major axis. We see, for highly inclined cases (is > 30°), the 
evolution of planetesimal semi-major axis is nearly independent of the eccentricity and nodal precession of the gas. We thus feel justified 
in using the simple CGD model for the remainder of our investigation due to the fact that we are concentrating on the high inclination 
cases. A detailed investigation of the planetesimal dynamics, including the reason for large inward jumps is provided in §|3] 
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Fig. 4. — Schematic diagram of different regimes for the damping of planetesimal semimajor axis. The damping is dominated by r) 
(the sub-Keplcrian gas velocity) in region S, the evolution of i and e dominates the gas damping in region B. Region B is divided into 3 
sub-regions depending on the particular mechanism by which e and i are excited; secular perturbations (Bl), Kozai effect (B2), and other 
effects such as mean motion resonance by the companion star (B3). SB is a transition region between S and B, in which contributions from 
e, i and r\ can all be significant. The dashed line divides the plot into two parts: e has more contribution to the damping of sem i-major 
axis than than i does in the top-left part, while in the bottom-right part the opposite is true. For full details, see the text in S I3.1I 
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Fig. 5. — The evolution of orbital eccentricity (e, dashed, top), inclination (i, solid, top), velocity relative to the gas (V re i in units of local 
Kcplcrian velocity , dashed, middle), gas drag force (F$, solid middle), semimajor axis (red solid, bottom) of a planctesimal with a p o = 7 
AU, and R p = 5 km under the coupled effects of perturbations from a binary (as = 50 AU, es = 0.3, is = 30°, Ma = 2Mb = Mq) and 
gas drag from a gas disk (1 MMSN) . Planetesimal orbital eccentricities and inclinations are initially set as very small values (< 10 — 4 ), 
and their other initial angular orbital elements are set as random values from to 360° . In the bottom panel, we also plot the results for 
planctcsimals with different initial semi- major axis a p o- The dashed curves are the results of retrograde case (is = 150°) for comparison. 
We see that all planetesimals' inward migrations transit from a fast regime to a slow o ne a t a similar turn-over semi-major axis, which is 
independent of a p o- This turn-over semi- major axis can be well described with a c (Eg 1 171 and see also the horizonal dot-dash line in the 

bottom panel). Note that i is measured relative to the plane of the gas disk rather than to the binary orbital plane. The vertical dot-dash 
line denotes the time at the turnover point for the case of a p o = 7 AU. 
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Fig. 6. — The same as Fig[5] but for the Kozai-on case with ig = 60°. We see that high i and e (induced by the Kozai effect) significantly 
increase the inward migration; planctcsimals can even "jump" into the very inner region, and unlike to that in the Kozai-off case, we now 
find that the turn-over semi-major axis (where the jumping stops) is no longer independent on a p o (see FigO. The horizonal dot-dash 
line denote i = ig = 60° . 
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a R = 50 AU, e R = 0.3, i R =60°, M A =2M R = Msu n 




/ 



/ o 




2 4 6 8 10 

a p0 (AU) 

Fig. 7. — The turn-over semi-major axis for the high inclination Kozai-on case (denoted with "y" here) v.s. a p o for the case with as = 50 
AU, &b = 0-3, «s = 60°, and Ma = 2Mg = Mq. Planctcsimals' orbital eccentricities and inclinations are initially set as very small values 
(< 10~ 4 ), and their other initial angular orbital elements are set as random values from to 360°. We see that y ~ a p o if a p o < a c , i.e., 
a p0 IS 1-5 AU. For 1.5 < a p o < 8 AU, y follows an oscillatory pattern, whose amplitude varies between a lower limit of a c (; om ) ~ 0.3 AU 
and an upper limit of a c , i.e., y ~ 0.30 — 1.5 AU. For a p o beyond ~ 8 AU, the results become chaotic, but the upper and lower limits 
(y ~ 0.30 - 1.5 AU) still hold. 
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Fig. 8. — Empirical fitting of a c n ow \ as a function of as for cases of %b = 50° (dot, open circle), 60° (dashed-dot, open triangle), and 
70° (dashed-dot-dot, open upside down triangle). The circle, triangle symbols are the results of individual measurements made by plotting 
fig ures similar to Fig[T] We also plot the disk stability boundary (solid line), a d ~ 0.22ag for Ma = 2Mg = ^© an d e B = 0.3 according 
to lHolman fc Wic gcrt (f99!j) , and a c (dashed) as references. We find that a good empirical fits for a c (; „) can be given as Eq |18l In the 
inset figure at the bottom right, wc also plot the measured a, a n ow \ as a function of ig. The open symbols are for prograde cases while 
the filled ones are for the corresponding retrograde cases, while the i = 90° case is plotted using an asterisk. We see that the results are 
symmetrical about i = 90° . 
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Fig. 9.— Planetesimal evolution for the cases in which M A = 2M B = M Q , a B = 50 AU, e B = 0.3, i B = 10° , 30° , 50° , 60° , 70° , 80° , 
R p = 5 km The second panel in the right-hand column (highlighted with blue bold axes) shows the results of the standard case (§ 14.21 1. We 
find in highly inclined cases that the planetesimal disk surface density is highly enhanced (£/£rj ~ 10) in the inner region within a c , while 
i t is sig nificantly depicted beyond a c . In addition, all results seems to be approximately symmetrical about i = 90°. For details see text in 
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Fig. 11.— Planetesimal evolution for the cases in which M A = M Q , q B = M B /M A = 0.2,0.5,1.0, a B = 50 AU, e B = 0,0.3,0.6, 
i B = 60°, R p = 5, 100 km. The middle panel (hi ghligh ted with blue bold axes) shows the results of the standard case (§ 14,21 1, The specific 
effects of q B , e B , R p and fd are summarized in § 14,3,31 In all cases, the general result holds, i.e., disk surface density is highly enhanced in 
the inner region within a c , while it is significantly depleted beyond a c . 
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Fig. 12. — Analytical estimate (Eq[2TJ of surface density enhancement (£/£) as a function of semimajor axis of binary systems with 
fixed orbital eccentricity eg = 0.3 and mass ratio Mb/Ma = 0.5. Results of four cases with different initial disk profiles are plotted in lines 
of different styles. The two vertical dash-dots lines, ag ~ 12.3 AU and ag ~ 135.8 AU mark the turn-over points where = dice = 2.7 
AU and a c = a; ce = 2.7 AU, respectively. 
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Fig. 13. — Evolution of relative velocities (top) and semi-major axis (bottom) for 10,000 planctcsimals in a system with the standard 
binary configuration, Ma = Mq, qs = Mb /Ma = 0.5, ag = 50 AU, eg = 0.3, %b = 60°. Planctcsimals with a centered Gaussian size 
distribution between R p = 1 km and R p = 10 km, are initially randomly distributed between 0.3 AU and 10 AU (flat distribution) with 
near circular and coplanar orbits. In the top panel, the contours denote the rate of collisions (in arbitrary units), and the thick solid line 
denotes the evolution of average relative velocity, for all collisions we have searched between 0.3 AU and 10 AU. In the bottom panel, the 
contours denote the changes in the planetesimal surface density, i.e., E/£o- We see that with the planetesimals piling up in the inner region, 
their relative velocities and collision rates are significantly reduced and increased respectively, which thus favors planetesimal growth in the 
inner region. 
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APPENDIX 

MAXIMUM FORCED ECCENTRICITY 

If the secular perturbation (Kozai effect off) dominates the orbital evolution of a planetesimal. Then the average 
eccentricity of the planetesimal can be estimated as its forced eccentricity e/, where 

4 as 

which is proportion to the plane tesimal semima j or ax is. Thus the maximum ef is obtained at the edge of the 
circumprimary disk. According to iPichardo et aTl ()2005h . the radial size of the circumprimary disk can be estimated 
as 

Rd ~ i?£ 99 0.733(l - es) 1 - 20 ^- 07 (A2) 
, where RE gg is the radius of the Roche lobe calculated by Eggleton (1983), 

0.49 gi 

R Egg ~ 2/3 1/3 a B . (A3) 

0.6^' +m{l + q 1 / ) 

Here, qi = Ma/M b and g2 = Mb / {Ma + Mb) are two mass ratios, and Ma and Mb are the masses of the primary 
and the secondary respectively. For a typical binary mass ratio, gi = 1/2, we have Rd ~ 0.38as(l — es) 1 ' 2 ■ Substituting 
this Rd into Ea lAll then the maximum forced eccentricity is e/ ~ 0.1. 



